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I.   GENERALIZED  LEAST  SQUARES 

Consider  the  linear  model: 
y  =  X3  +  e  , 

2 

where  e  is  N(0,cr  z)   and  X  is  an  txp  matrix  of  rank  p. 

The  ordinary  least  squares  (OLS)  estimate  of  0  for  this 
model  is: 


-1 
0  =  (X'X)   X'y. 


Furthermore,  an  estimate  of  Cov(0)  is  given  by; 


^   ^         -12 
Cov(0)  =  (X'X)    S  , 


2  -1  -1 

where  S   =(n-p)   y'(I-X(X'X)   X')y. 


OLS  is  appropriate  for  the  linear  model  presented 

above  if  the  covariance  matrix  of  the  error  term  (e) 

conforms  to  particular  patterns.   The  necessary 

conditions  are  discussed  in  a  paper  published  by  Huynh 

and  Feldt  (1970).   According  to  Milliken  and  Johnson 

(1984),  the  Huynh-Feldt  (H-F)  condition  is  equivalent  to 

specifying  that  the  variances  of  the  differences  between 

pairs  of  errors,  such  as  e   -  e    ,  are  equal  for  all 

i     i' 


i  and  i',  i  ^  1'.   If  a  model  has  more  than  one  error 
term,  the  covariance  matrices  of  each  error  vector  must 
satisfy  the  H-F  condition  for  OLS  to  be  valid. 

There  are  essentially  3  consequences  associated 
with  using  OLS  when  the  H-F  condition  does  not  hold. 
They  are: 

1.   (3  is  unbiased.   E(0)  =  3. 


2.  Cov(3)  will  be  biased. 

3.  Ordinary  least-squares  will  give  inefficient 
predictions. 


Departures  from  the  H-F  condition  can  be  remedied 
using  generalized  least  squares  (GLS).   The  GLS  estimate 
of  S  is: 

-1   -1    -1 
3  =  (X'z:   X)   X'E   y, 

which  can  be  computed  as  a  regression  of  Ay  on  AX,  where 

-1     • 
E    =  A  A  . 

An  estimate  of  Cov(3)  is  given  by: 


A     ^ 


-1   -1   2 

cov(3)  =  (x'e:  X)   s  , 

2  -1         -1         _1  _1     _1       _;l 

where  S   =  (n  -  p)   y' (e    -  e   x(X'e   X)  X'e      )y. 
According  to  Arnold  (1981),  the  following  are 


properties  of  the  GLS  estimators  }i    (where  ^  =  X|3)  and 
2 


"      2 
1.   (M,S  )  is  a  complete  sufficient  statistic 


2 

2.  p.   and  S   are  the  minimum  variance  unbiased 

2 
estimators  of  }i   and  a  . 


2 
3.   fi  and  ((n  -  p)/n]S  are  the  maximum  likelihood 

2 

estimators  of  }i   and  a  . 


SIMPLE  LINEAR  REGRESSION  MODEL 
WITH  AUTOREGRESSIVE  ERRORS 


To  illustrate  the  problems  that  arise  when  serially 
correlated  disturbances  exist,  let  us  postulate  the 
following  model: 


Y=0+0x+e,    k  =  1, ,t, 

k    o    Ik    k 


where  we  assume  the  error  term  e   to  be  first-order 

k 


autoregressive  and  of  the  for 


e  =  *e      +  e 
k     k  -  1    k 


m: 


where  *  represents  the  coefficient  of  correlation  and 


hence  -1  <  *  <  1.   Furthermore,  e  satisfies  the 

k 
following  assumptions  for  all  k: 


E(e  )  =  0 
k 


E(e  e      )  =  a 
k  k  +  r      e 

=  0 


r  =  0 


r  ^  0 


Given  the  above  assumption's,  the  var iance-covariance 
matrix  for  the  serially  correlated  disturbances  can  be 
written  as: 


2      2 
E(eE  •  )  =  a  V  =  a 
e      e 


•» 


* 


* 


« 


t-1    t-2    t-3 
♦      ♦      ♦ 


t-1  , 


•» 


t-2 


.  * 


.  1 


txt 


2     2  2 

where  a   =  a   /  (1  -  *  )   for  all  t 
e    e 


CONSEQUENCES  OF  AUTOCORRELATED  ERRORS 

In  his  book  Econometric  Methods  (2nd  ed .  )  ,  Johnston 
describes  three  consequences  of  applying  ordinary  least- 
squares  to  a  model  with  disturbances  that  are 
autocorrelated.   They  are: 


1.  the  estimates  of  |3  are  unbiased, 

2.  the  sampling  variances  of  the  regression 
coefficients  are  likely  to  be  seriously 
underestimated,  and 

3.  ordinary  least-squares  will  give  inefficient 
predictions. 


Because  of  the  aforementioned  concerns,  any  t  or  F 
test  performed  on  the  parameters  of  the  model  will  be 

2 
incorrect  because  cov(0)  is  biased  and  because  SSE/a   is 

e 
not  distributed  as  a  chi-square  nor  is  it  independent  of 
(b  -  0),  (where  b  is  the  OLS  estimate  of  0). 

Given  that  the  matrix  L    is  known,  the  generalized 
least  squares  estimator  of  0  and  its  corresponding 
standard  error  can  be  computed  so  that  the  usual 
confidence  interval  and  hypothesis  test  of  0  can  be 
constructed.   The  minimum  variance  unbiased  estimator  of 
0  in  this  situation  is  provided  by  generalized  least- 
squares  (GLS). 


GLS  ESTIMATES 

Consider  the  matrix  representation  of  the  simple 
linear  regression  model: 

y  =  X0  +  e 


2 

with  E(e)  =  0,  E(ee' )  =  a  c. 

e 


Given  that  the  error  term  is  known  to  be 
autoregressive  and  that  *  is  known,  GLS  estimates  for 
this  model  can  be  obtained  either  directly  or  by  a  two- 
stage  procedure. 

The  direct  GLS  estimator  of  (3  is: 

-1   -1   -1 
0  =  (X'E   X  )   X'E   y. 

The  variance  of  this  estimator  is: 


"   "     2    -1   -1 
VAR(0)  =  3  (X'E   X)    , 

2-1 
where  s   =  (J'e  d/(n  -  2)   and   <5  =  y  -  X3. 

The  two-stage  procedure  consists  of: 

1.  transforming  the  elements  of  the  Y  vector  and 
the  X  matrix  by  appropriate  functions  of  *,  and 

2.  applying  OLS  to  the  transformed  data. 


The  model  Is  transformed  by  premultlplying  the  Y  vector, 
X  matrix  and  error  vector  by  a  matrix  A  (Ay  =  AX0  +  Ae) 

2 

such  that  E(Aee'A')  =  a  I.   The  appropriate  transform- 

e 

atlon  matrix  in  the  context  of  a  simple  linear 
regression  model  with  t  observations  and  an  AR(1)  error 
structure  is: 


A  = 


SQRTd  -  *  )    0    0 

-  ♦         10 

0      -  *    1 


0  0  0 
0  0  0 
0    0    0 


0 
0 


0 
0 


0  .  . 
0  .  . 


.  -  *    1    0 
0  -  *    1 


txt 


Ordinary  least-squares  is  applied  to  the  following 
transformed  data: 


Ay  = 


1 


SQRTd  -  *  )Y 


Y  -    *Y 

2  1 

Y  -  ♦Y 

3  2 


Y   -  «Y 

t      t  -  1 


AX  = 


2  2 

SQRTd  -  *  )    SQRTd  -  *  )X 


1  -  * 


1  -  * 


X   -  *X 
2      1 


X   -  *X 

t      t  -  1 


The  same  transformation  applied  to  the  simple 
linear  regression  can  be  extended  to  a  repeated  measures 
experimental  design  for  which  the  subplot  errors  are 
serially  correlated. 


II.   REPEATED  MEASURES  WITH  AUTOREGRESSIVE 
SUBPLOT  ERRORS 


A  repeated  measures  model  is  similar  to  a  split- 
plot  model,  except  for  the  assumptions  inherent  to  the 
respective  models.   The  model  is: 


'=M+a   +it    +T   +  (aT)    +  e 
^^^       i     ij    k       ik    ijk 


i  =  l,...,a,  j  =  l,...,n,  k  =  l,...,t. 
Where  }i    is  the  overall  mean. 


(1), 


■^  represents  whole  plot  treatment  effects 


n       Is  the  error  term  associated  with  the 
Ij      whole  plots. 


(aT)      represents  the  whole  plot  treatment  by 
ik   subplot  treatment  interaction. 


T        represents  subplot  treatment  effects 
j       (usually  time  for  a  repeated  measures 
experimental  design). 


e        is  the  error  term  associated  with  the 
Ijk     subplots. 


The  form  of  the  covariance  matrices  associated  with 
the  error  terms  determines  the  appropriate  analysis  of  a 
repeated  measures  experimental  design.   Two  conditions 
will  be  discussed  in  particular: 


1.  the  case  for  which  a  split-plot  analysis  is 
appropriate,  and 

2.  the  case  where  the  covariance  matrix  of  the 
subplot  errors  is  assumed  to  be  autoregressive 


SPLIT-PLOT  ANALYSIS 

A  split-plot  analysis  of  a  repeated  measures 
experimental  design  is  appropriate  if  the  covariance 
matrix  for  each  error  term  in  the  model  conforms  to  the 
Huynh-Feldt  (H-F)  condition.   An  Analysis  of  Variance 
table  can  be  written  as: 


10 


Table_l ANOVA 

Source  of  Variation df SS 

«  a  -  1  SS(a) 

Error(Whole  Plot)       a(n  -  1)         SSE(Whole  Plot) 
T  t  -  1  SS(T) 

a  *  T  (a  -  l)(t  -  1)    SS(a  *  T) 

Error(Subplot)  a(n  -  l)(t  -  1)   SSE(Subplot) 


Source  of  Variation MS. 


ANOVA 


SS(a)/(a  -  1)         MS(a) 

MSE(Wole  Plot) 


Error(Whole  Plot) 


Error (Subplot ) 


SS(T)/(t  -  1)         MS(T) 

MSE(Subplot) 


a  *  T  SS(a  *  T)  HSioL    *    T) 

(a  -l)(t  -1)       MSE(Subplot) 


AUTOREGRESSIVE  REPEATED  MEASURES 

If  the  subplot  error  structure  is  autoregressive  by 
nature,  the  test  statistics  computed  from  a  split-plot 
analysis  of  the  subplot  treatment  (T)  and  the  whole  plot 
treatment  by  subplot  treatment  interaction  (a  *  T)  is 


11 


erroneous.   However,  the  whole  plot  tests  still  remain 

valid  even  if  the  subplot  error  (e    )  is  not 

ijk 

independent. 

Three  methods  are  proposed  to  remedy  the 
shortcomings  of  ordinary  least  squares  (split-plot 
analysis)  when  the  subplot  error  structure  is 
autoregressive.   They  are: 


1.  Box's  adjusted  degrees  of  freedom  associated 
with  the  subplot  analysis, 

2.  generalized  least  squares  (GLS),  and 

3.  generalized  linear  model  -  generalized  least 
squares  (GLM  -  GLS) . 


Following  the  introduction  of  each  methodology,  an 
example  will  be  presented  to  illustrate  each  concept. 

BOX'S  ADJUSTED  DEGREES  OF  FREEDOM 

According  to  Box  [1954],  an  adjustment  can  be  made 

to  the  degrees  of  freedom  of  an  experimental  design  with 

2 
an  error  structure  of  the  form  a  E.   The  degrees  of 

freedom  are  multiplied  by  a  constant  which  is  computed 
as  follows: 


2  -      -    2 
t  (V    -  V   ) 

aa ,_^,__ 


ab     2      a-2     2- 
(t  -  1)  (E  E  V     -  2t  E  V     +  t   V 
1  1  ab       la. 


12 


where  (V   )  =  V  is  the  covarlance  matrix  o£  the  t 
ab 

observations  on  a  subject.   V    is  the  mean  of  the 

aa 

diagonal  elements  of  V,  V    is  the  mean  of  all  elements 


of  V  and  V   is  the  mean  of  the  elements  in  row  "a"  of 
a. 

of  V. 


Wallenstein  and  Fleiss  [1979]  have  extended  Box's 
idea  to  encompass  situations  where  the  subplot  errors 
are  autoregressive  by  nature.   For  such  a  covarlance 
structure,  the  degrees  of  freedom  associated  with  the 
subplot  treatment,  the  Interaction  between  the  whole 
plot  treatment  and  the  subplot  treatment,  and  the 
subplot  error  are  multiplied  by  the  constant: 

Y      =     R/B,      where 


2  t      2 

R    =    (1    +   *)    {t(l    -   *)[t(l    -*)-(!+   *)]    +    2*(1    -   *    )} 

and, 

t  t 

B    =    (t    -    1){4*(1    -*)(!    +    *)[*{!    +   *)(i    -   4,    ) 

2  t 

+   t(l   -   *)(i    +   «)      -   t*(l    -   *)(i    +   *    )] 

2  22  2t  4  2t+l 

-    t    (1    -   ♦)     [2#    (1    -   *      )+(!+*)       +    4(1    +    *)    *  ] 

3  3  2 

+   t    ( 1    -   ♦ )    (1    +   * )  ( 1    +   ,^    ) } 
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The  degrees  of  freedom  to  be  adjusted  can  be 
obtained  from  the  split-plot  ANOVA  table  above.   The 
appropriate  adjustments  are  presented  in  Table  2  below. 


Table  -J,        Actual  and  Adjusted  Degrees  of  Freed 


om 


Source_of_Variation ^Actual_df ^Adjusted  df 


* 


a 


a  -  1  (a  -  1) 


Error(Whole  Plot)       a(n  -  1)  a(n  -  1) 

'''                      t  -  1  Y(t  -  1) 

^  *  "^                 (a  -  l)(t  -  1)  Y[(a  -  l)(t  -1)1 

Error(Subplot)  a(n  -  l)(t  -  1)  Y[a(n  -  l)(t  -1)] 


*  Unchanged 

The  appropriate  F  -  statistics  can  be  computed  once 
the  adjusted  degrees  of  freedom  have  been  computed  for 
the  subplot.   For  example,  assuming  the  subplot 
treatment  is  time  (T),  the  appropriate  computation  is: 

F  =  MS(Time)/MS(Error(Time)) 

where  MS(Time)  =  SS(Tlme)/Ad j .  df  Time,  and 
MS(Error(Time))  =  SS(Error (Time) )/Adj .  df  Error(Time). 

Since  the  numerator  and  denominator  sums  of  squares 
have  been  divided  by  degrees  of  freedom  that  have  been 
adjusted  by  the  same  factor, y,  the  F-statistic  for  the 
time  effect  will  be  exactly  the  same  as  the  F-statistic 
given  in  the  standard  split-plot  analysis  above. 
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However,  the  critical  region  used  In  the  analysis 
will  differ  from  that  of  the  standard  split-plot 
analysis. 

Table  3  below  gives  values  of  Y  for  *  =   0.05, 
0.10,  0.15,  0.20,  0.25,  0.30,  0.35,  0.40,  0.50,  0.55, 
0.60,  0.65,  0.70,  0.75,  0.80,  0.85,  0.90,  0.95,  and  for 
t  =  3,  4,  5,  7,  10,  25,  50,  and  infinity. 

Table  -_3;  values  for  gamma  when  the  correlation  matrix 

HAS    A    SIMPLEX    PATTERN 
NUMBER    OF    TIME   POINTS 
♦  3  4  5  7  10  25  50  <x> 

0.05  0.999    0.998    0.997    0.997    0.996    0.995    0.995    0.995 

0.10  0.996    0.992    0.990    0.987    0.985    0.982    0.981    0.980 

0.15  0.991    0.983    0.977    0.971    0.966    0.960    0.958    0.956 

0.20  0.985    0.971    0.961    0.949    0.940    0.929    0.926    0.923 

0.25  0.977    0.956    0.941    0.923    0.909    0.892    0.887    0.882 

0.30  0.968    0.940    0.919    0.892    0.872    0.848    0.841    0.835 

0.35  0.958    0.921    0.893    0.858    0.831    0.799    0.790    0.782 

0.40  0.948    0.901    0.866    0.822    0.787    0.745    0.734    0.724 

0.45  0.936    0.880    0.838    0.783    0.740    0.688    0.674    0.663 

0.50  0.925    0.859    0.809    0.743    0.692    0.692    0.613    0.600 

0.55  0.912    0.837    0.779    0.703    0.643    0.568    0.550    0.536 

0.60  0.900    0.814    0.749    0.663    0.594    0.507    0.486    0.471 

0.65  0.887    0.791    0.719    0.623    0.546    0.446    0.423    0.406 

0.70  0.875    0.769    0.689    0.584    0.449    0.387    0.360    0.342 

0.75  0.862    0.747    0.660    0.546    0.454    0.330    0.299    0.280 
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VALUES  FOR  GAMMA  WHEN  THE  CORRELATION  MATRIX 
HAS  A  SIMPLEX  PATTERN 

NUMBER  OF  TIME  POINTS 

♦     3      4      5      7     10     25     50      « 

0.80      0.849    0.724    0.631    0.510    0.411    0.275    0.240    0.220 

0.85      0.837    0.703    0.604    0.475    0.371    0.225    0.185    0.161 

0.90      0.824    0.682    0.577    0.442    0.333    0.179    0.133    0.105 

0.95      0.812    0.661    0.551    0.410    0.298    0.138    0.088    0.051 

LOWER 

BOUND  0.800  0.641  0.526  0.381  0.266  0.103  0.051 


(T-1)    0.500    0.333    0.250    0.167    0.111    0.042    0.020 


Source:  REPEATED  MEASUREMENTS  ANALYSIS  OF  VARIANCE  WHEN 

THE  CORRELATIONS  HAVE  A  CERTAIN  PATTERN  By 

Sylan  Wallenstein  &   Joseph  1.  Fleiss  (1979:  p.   231) 


Table  3  is  useful  in  determining  whether  an 
adjustment  is  needed  to  the  subplot  analysis.   Box's 
recommended  adjustment  to  the  degrees  of  freedom  does 
not  result  in  tests  that  are  as  powerful  as  tests 
generated  by  generalized  least  squares  (GLS).   However, 
GLS  is  time  consuming.   If  the  adjustment  factor  does 
not  change  the  degrees  of  freedom  much,  one  need  not 
expend  the  time  nor  effort  performing  a  GLS  analysis  on 
a  repeated  measures  experiment  with  an  autoregressive 
error  structure. 

GLS  ESTIMATES 
An  alternative  expression  for  the  split-plot  model 
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which  emphasizes  the  nature  of  the  whole  plots  is 


Y     =   p.         +  T   +  (aT)    +  e      (2), 
ijk    ij    k       ik    ijk 

i  =  l,...,a,     j  =  l,...,b,    k  =  l,...,t. 


where  p.  is  a  block  effect  for  the  whole  plot 

ij 

experimental  units  and  is  equal  Ji  +  a  +   n      .      The 

i     ij 

blocks  average  out  the  effect  of  the  autoregressive 

subplot  disturbances,  validating  the  whole  plot  analysis 

of  a  split-plot  model. 

Model  two  above  can  be  represented  in  matrix 
notation  as: 

y  =  X0  +  e  , 


where  y  is  an  (ant)xl  response  vector  whose  members 
consist  of  y   , 
ijk 


X  is  an  (ant)x[a(n+t)  +  t)]  design  matrix,  and 

0  is  an  [a(n+t)  +  t]xl  vector  of  parameter 
estimates.   3'  =  t^',  T',  aT']  , 


where   W    =    [}i     ,    n      ,    .  .  .    ,    p.      ]•  , 

11   12         an 


T'  =  [T  ,  T  ,  . . .  ,  T  ]•  ,  and 
12         t 


aT'  =  [aT   ,  aT   ,  . . .  ,  aT   1  ' 
11    12  at 
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E  Is  an  (ant)xl  error  vector  whose  members  consist 
of  e 

ijk 


To  illustrate  what  the  design  matrix  for  such  a 
model  looks  like,  let  us  look  at  the  simplest  case 
(a=n=t=2).   Assume  an  experiment  was  conducted  to 
monitor  the  effect  of  a  drug  on  the  heart  rate  of  a 
patient.   The  experiment  was  designed  so  that  two 
patients  (n=2)  were  randomly  assigned  to  each  of  the  two 
drugs  (a=2)  tested.   The  heart  rates  per  minute  were 
obtained  for  each  of  two  weeks  (t=2).   The  design  matrix 
for  this  scenario  is: 


X  = 


1 

0    0    0 

1 

0 

1 

0 

0 

0 

1 

0    0    0 

0 

1 

0 

1 

0 

0 

0 

10    0 

1 

0 

1 

0 

0 

0 

0 

10    0 

0 

1 

0 

1 

0 

0 

0 

0    10 

1 

0 

0 

0 

1 

0 

0 

0    10 

0 

1 

0 

0 

0 

1 

0 

0    0    1 

1 

0 

0 

0 

1 

0 

0 

0    0    1 

0 

1 

0 

0 

0 

1 

where  the  first  partition  consists  of  indicator 
variables  representing  the  whole  plots,  the  second 
partition  represents  time,  and  the  final  partition  of 
indictor  variables  represents  the  drug  by  time 
interaction. 


Since  the  design  matrix  for  the  split-plot  model  is 
not  full  column  rank,  a  generalized  inverse  needs  to  be 
used  to  estimate  0.   For  example,  given  that  e  is 
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2 

N(0,a  I): 

E 


0  =  (X'X)  X'y 

Given  that  e  is  autoregressive,  the  covariance  matrix 
for  the  subplot  error  is: 


2      2 
E(E'E)   =  a   E  =  a   (I  ®  V) 
E      £   an 


The  following  steps  should  be  taken  to  achieve  an 
appropriate  analysis  of  a  repeated  measures  experimental 
design  with  autoregresslve  subplot  disturbances: 

1.  Perform  a  split-plot  analysis  using  OLS.   The 
F-statistics  for  the  whole-plot  analysis  will 
be  correct,  however,  the  subplot  analysis  will 
give  erroneous  F-statistics.   The  next  two 
steps  are  needed  to  correct  the  subplot 
analysis . 

2.  The  second  step  consists  of  premultiplying 
the  design  matrix,  X,  and  the  response  vector, 

y,  by  the  matrix  Z,  where  Z  =  I   ®A  . 

an 


2 
3.   Because  var(ZE)  =  cr   I  ,  one  can  obtain  GLS 

E 

2 
estimators  of  (3  and  a  by  running  OLS  on  the 
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variables  o£  the  model: 

*    *      * 

y  =  X  0  +  £ 

* 
where  y  =  Zy, 

* 
X  =  zx, 

* 

e  =  ze 

The  output  from  regressing  Zy  on  ZX  will 
provide  appropriate  F-statistics  for  the 
subplot  part  of  the  analysis. 

The  OLS  estimator  of  (3  (GLM  estimator  of  (3  for 
the  split-plot  model)  is: 

•  -1   -   •  -1       II     _    I  I 
0   =(XE   X)   XE   y=(XZZX)    XZZy 

Recall  that  GLS  produces  minimum  variance  unbiased 
estimators  of  the  mean  vector  and  covariance  matrix. 
However, the  design  matrix,  X,  needs  to  be  constructed  to 
implement  this  procedure.   This  is  time  consuming. 

GLM-GLS 

The  procedures  and  properties  for  the  generalize 
linear  model  -  general  least  squares  (GLM-GLS)  method  is 
similar  to  that  of  the  regular  GLS  analysis  presented 
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above.   However,  GLM-GLS  is  less  time  consuming, 
Consider  the  transformation: 


+ 
y     =(l-*)y      ,        k=l 
ijk  ijk 


ijk     ij(k-l) 


In  terms  of  the  split-plot  model,  this  translates  into: 


+      +     +     +     + 
y    =  ]x        +T+aT   +e     (3), 
ijk    ij    k     ik    ijk 


+ 
where   )k        =  (1  -  *)^    ,      k  =  1,  ...,  t 
ij  ij 


+ 
T=(1-*)T    ,       k=l 
k  k 

=  T  -  ♦?     ,      k  =  2,  ...,  t. 
k     k-1 

+ 
aT    =  (1  -  *)  aT        ,    k  =  1 
ik  ik 

=  aT    -  *aT        ,    k  =  2,  ...,  t, 
ik       i(k-l) 

+ 

e    =(l-»)e         ,k=l 
ijk  ijk 

ijk 


The  errors  from  this  transformed  model  are 


now 
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Independent,  but  do  not  have  constant  variance 


Var(E    ) 
ijk 


(1  -  i) 


1  -  * 


k  =  1 


cr         ,   k  =  2, . .  .,  t. 

e 

However,  a  nonconstant  variance  can  easily  be  dealt  with 

by  most  computer  packages. 

Model  three  can  be  represented  in  matrix  notation  as: 


+     +    + 
y  =  X|3   +  E 


where. 


+ 
0   = 


aT 


+ 
M  = 


(1  -  *)  M' 


DT' 


aT   =  (I  X  D)  aT' 

a 


where. 


D  = 


(1  -  *)   0    0 

-  *      1    0 

0    -  *    1 


0     0     0  ,  . 
0   '  0     0  .  . 


The  GLS  estimator  of  3   is  given  by: 


^+  -    + 

3   =  (X'WX)  X'Wy 


where  W  =  I   ®  G  ,  and 
an 


-  * 
0 
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0  0  0 
0  0  0 
0     0    0 


1    0 

*   1 


txt 


G  = 


2 

1  - 

-   * 
2 

0 

0       .     . 

0 

0 

0 

(1  - 

-   ♦) 

0 

1 

0       .     . 

0 

0 

0 

0 

0 

1    .   . 

0 

0 

0 

0 
0 


0     0 
0     0 


0    10 
0    0    1 


txt 


Now,  D  is  a  non-singular  invertible  matrix  with  inverse: 
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-1 


(1 

t 

*) 

0 

(1 

♦  ) 

1 

JL 

> 

« 

•    •    • 


(1  -  ♦) 


t-1     t-2    t-3    t-4 
J6 *        *        *  » 


(1  -  *) 


txt 


so  that  a  test  of   H  :  T  =  0   is  the  same  as  testing 

o 


H  :  T  =0,  and  a  test  of  H  :  aT  =  0   is  the  same  as 
o  o 


testing   H  :  aT   =0.   As  with  GLS,  the  whole-plot 
o 

F-tests  produced  by  a  standard  split-plot  analysis  will 
be  correct. 

A  statistical  package  such  as  SAS  will  construct 
the  design  matrix,  X,  and  the  matrix  of  weights,  W. 
One  need  only  specify  the  weights  along  the  diagonal  of 


the  matrix  W  and  construct  y  ,  which  is  relatively 
easy. 
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Presented  below  is  an  example  of  the  SAS  code 
needed  to  yield  an  appropriate  analysis  of  the  subplots 
for  the  simplest  case  presented  above  (a=n=t=2): 


DATA  EXAMPLES- 
INPUT  PERSON  TIME  DRUG  Y; 

/*  COEFFICIENT  OF  CORRELATION  IS  ASSUMED*/; 
/*  FOR  ILLUSTRATIVE  PURPOSES  TO  EQUAL  0.6  */; 
RHO  =0.6; 

/*  THE  WEIGHTS  FOR  THE  MATRIX  W  NEED  TO  BE  */; 
/*  SPECIFIED  AND  ARE  ASSUMED  TO  BE  4  FOR  */; 
/*  TIME  PERIOD  1  AND  1  FOR  TIME  PERIOD  2  */; 
IF  TIME  =  1  THEN 

W=(l  -  RH0"2)/(1  -  RH0)"2; 
ELSE  W=l; 
LAGY  =  LAG(Y); 
IF  TIME  =  1  THEN 

YPLUS  =  (1  -  RHO)*Y; 
ELSE  YPLUS  =  Y  -  RHO*LAGY; 
CARDS; 


data 


PROC  GLM  DATA=EX AMPLE; 
CLASSES  PERSON  TIME  DRUG; 

MODEL  YPLUS  =  PERSON (DRUG)  TIME  DRUG*TIME/ 
NOINT  SOLUTION  SSI; 
WEIGHT  W; 


CONTRASTS 


The  GLM-GLS  parameter  estimates  of  the  subplot 
treatment,  usually  time,  are  linear  functions  of  the 
GLM  parameter  estimates  of  the  subplot  treatment.   The 
relationship  Is  given  by: 
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+ 
T   =  DT. 

+ 
Conversely,  T  is  a  function  of  T  .   The  relationship  is 


given  by: 


-1  + 
T  =  D  T  . 


Given  the  relationship  that  exists  between  T  and 

+ 
T  ,  any  contrasts  of  the  form  CT,  where  C  is  a  matrix  of 

desired  contrasts,  can  be  estimated  alternatively  as 

-1  + 
CD  T  . 


RECOMMENDED  STRATEGY 


The  following  strategy  is  recommmended  for  experimental 
designs  with  autoregressive  subplot  errors: 

1.  Fit  the  model  by  OLS .   From  this  analysis: 

-  test  the  whole  plot, 

-  estimate  ♦  with  the  residuals  obtained 
from  OLS  (Albohali,  1983),  and 

-  consult  Table  3  above  to  see  if  a  GLS  or 
GLM-GLS  analysis  needs  to  be  performed 
on  the  subplot.   If  the  adjustment 
factor,  Y,  is  close  to  one,  a  standard 
split-plot  analysis  will  be  appropriate 
for  both  the  whole  plot  and  the  subplot. 

2.  Do  the  necessary  transformations.  Keep  in 
mind  that  with  a  GLM-GLS  analysis  one  need 
not  produce  the  design  matrix, X. 
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III.   EXAMPLE 

BOX'S  ADJUSTED  DEGREES  OF  FREEDOM 

To  illustrate  Box's  adjustment  to  the  subplot 
degrees  of  freedom,  consider  an  experiment  in  which  the 
effect  of  a  dose  of  a  drug  on  the  growth  of  rats  is 
monitored.   Presented  below  in  Table  4  are  growth 
measurements  of  50  rats,  where  10  rats  were  randomly 
assigned  to  each  of  of  the  five  doses.   The  weights 
where  monitored  and  collected  for  11  weeks. 

Tabl^  -  ,4  Body  Weights  of  Rats 


WEEK 


DOSE   RAT    1    2    3    4    5    6    7    8    9    10    11 


10 


1  54  60  63  74  77   89   93   100  108  114  124 

2  69  75  81  90  97   120  114  119  126  138  143 

3  77  81  87  94  101  110  117  124  134  141  151 

4  64  69  77  83  88   96   104  109  120  123  131 

5  51  58  62  71  74   81   88   93   99  103  113 
64  71  77  89  90   100  106  114  122  134  139 

7  80  91  97  101  111  119  129  131  137  147  154 

8  79  85  89  99  104  105  116  121  132  139  147 

9  77  82  88  92  101  109  119  127  135  144  158 
79  84  91  98  107  114  119  131  137  146  155 


6 
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Week 


DOSE      RAT      123456789         10         11 


•5      1  62  71  75  79  87  91   100  105  111  121  124 

2  68  73  81   89  94  101  110  114  123  132  139 

3  94  102  109  110  128  133  147  151  153  171  184 

4  81  90  95  102  109  120  128  137  141  154  160 

5  64  69  72  76  84  89   97   103  108  114  124 

6  67  74  81  81  84  95   100  109  119  128  130 

7  73  80  86  89  97  101  110  116  117  135  141 

8  71  74  82  84  93  97   102  113  119  124  131 

9  69  74  79  89  94  100  107  113  124  134  139 
10  60  62  67  74  78  85   92   103  112  121  130 

1      1  59  63  66  75  80  87   99   104  110  115  124 

2  56  66  70  81  77  88   96   100  113  l20  130 

3  71  77  84  80  97  106  111  109  128  133  140 

4  59  64  69  76  85  88   96   104  110  119  126 

5  65  70  73  77  85  92   96   101  111  118  121 

6  61  69  77  81  89  92   107  111  118  127  132 

7  80  86  95  99  106  113  127  131  142  150  160 

8  74  80  84  90  99  101  108  117  126  133  140 

9  71  79  88  90  98  102  116  121  127  139  142 
10  69  75  80  86  96  97   104  113  122  129  138 

4      1  64  71  79  82  85  94   103  110  113  122  125 

2  53  57  61  72  74  76   81   91   99  100  105 

3  64  69  76  85  89  96   104  108  116  120  128 
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Week 
DOSE      RAT         123456789         10         11 


4  68  69  78  82  91  97   104  108  115  122  132 

5  69  74  80  85  90  99   104  114  123  129  133 

6  85  91  98  100  105  104  118  121  130  141  141 

7  75  82  85  92  99  107  112  125  130  137  146 

8  57  61  65  68  77  81   87   91   95  101  107 

9  69  72  77  80  84  91   96   103  109  116  125 
10  66  68  76  81  88  95   103  106  112  119  130 

1  57  64  70  76  80  90   93   99   105  113  118 

2  62  67  74  83  87  93   104  108  114  124  129 

3  60  68  73  80  83  94   101  106  112  122  131 

4  64  66  76  81  91  100  102  111  120  128  136 

5  57  60  67  73  67  64   75   85   89  98  105 

6  78  83  89  99  105  113  117  128  132  139  149 

7  81  81  92  100  108  119  120  133  138  149  157 

8  46  47  51  55  63  65   68   74   78  85  90 

9  69  72  74  76  77  82   82   90   95  101  103 
10  67  77  83  83  92  99   104  108  114  118  129 


Source:  Analysis  of  Messy  Data  Volumn  I:  Designed 
Experiments  by  Milliken  &   Johnson  (1984:  p.  371-372) 

The  appropriate  model  for  this  analysis  is: 

Y     =   }i   +  OL      +   K        +T+(aT)    +e 
i:Jk        I    ij    k       ik    ijk 

i  =  1,...,5,     j  =  l,...,io,      k  =  1,...,10. 
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Where  }i   Is  the  overall  mean. 

a      represents  dose  effect, 
1 

Tt      is  the  animal  error  term, 
ij 

T      represents  time  effect, 
J 


(aT)    represent  the  dose  by  time  interaction, 
ik 


e      is  the  time  interval  error  term  and  is 
ijk   considered  to  be  autoregressive . 


A  split-plot  analysis  of  the  data  yields  the 
following  Analysis  of  Variance  table: 


Tablg_? ANOVA 

Source  of  Var.     df gS MS 


Dose  4  10295.72  2573.93       1.53 

Error(Rat)  45  75668.30  1681.52 

Week  10  243381.13  24338.13    1783.51 

Dose  *  Week  40  1517.88  37.95      2.78 

Error(Week)  450  6140.80  13.65 


The  adjustment  in  this  case  is  made  to  the  degrees 
of  freedom  associated  with  week,  the  interaction  between 
dose  and  week,  and  the  subplot  error.   Each  is  multi- 
plied by  the  constant: 

Y  =  467.6421/811.1840  =  0.57649  (*  =  0.60). 
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The  adjusted  degrees  of  freedom  are  given  below; 


Table  -  6     Actual  and  Adjusted  Degrees  of  Freedo 


m 


Source  of  Var. 

Actual 

df 

Adiusted  d£ 

Dose 

4 

* 

4 

* 

Error(Rat) 

45 

45 

Week 

10 

5.7649 

Dose  *  Week 

40 

23.0596 

Error (Week) 

450 

259.4205 

*  Unchanged  (between  rat  comparison). 

The  F-statistics  in  the  split-plot  Anova  table 
should  be  compared  to  the  critical  regions  corresponding 
to  the  F  distribution  with  the  appropriately  adjusted 
degrees  of  freedom. 

The  appropriate  critical  regions  to  use  to  test 
for  a  significant  week  effect  and  to  test  for  a 
significant  dose  by  week  interaction,  respectively,  are: 


a  =  0.05 
F  =  2.15917 

(5.77,259.42) 


a  =  0.05 

F  =  1.57059 

(23.06,259.42) 


For  the  rat  data  above,  the  adjustment  factor  (v) 
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Is  substantially  below  an  adjustment  factor  of  one  (see 
table  3).   Hence  a  GLS  analysis  or  a  variant  thereof  is 
in  order.   The  magnitude  of  the  adjustment  factor  needed 
before  a  GLS  analysis  be  done  is  a  subjective  question. 


GLS 


To  illustrate  GLS,  let  us  take  another  look  at  the 
^   rat  data  presented  above.   The  first  step  involves 
running  a  split-plot  analysis  of  the  data  using  a 
procedure  such  as  SAS's  PROC  GLM.   The  F-statistic  from 
this  analysis  will  be  correct  for  the  whole  plot 
treatment,  dose  (see  split-plot  table  above). 

The  next  step  requires  the  creation  of  a  split-plot 
design  matrix  (X).  The  design  matrix  for  the  rat  data 
contains  121  variables,  five  for  dose,  fifty  for  whole 
plot  error  -  which  is  a  nested  error  structure 
consisting  of  rat  nested  within  dose,  eleven  for  week, 
and  fifty-five  for  the  dose  by  week  interaction.  A  no 
intercept  model  was  used. 

The  design  matrix  (X)  and  response  vector  (y)  is 

premultiplied  by  the  matrix  Z  =  I   ®  A,  where, 

an 
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A  = 


0.8  0  0 
0.6  1  0 
0    -  0.6   1 


0  0  0 
0  0  0 
0     0     0 


0      0     0  ...  -  0.6    1     0 
0      0     0  .  .  .    0    -  0.6   1 

11x11 
Note:  The  transformation  matrix,  Z,    is  a  550  x  550 
matrix. 

The  final  step  involves  running  OLS  on  the  121 
transformed  variables.   The  F-statistic  from  this 
analysis  will  be  correct  for  the  subplot  treatment, 
week,  and  for  the  interaction  term,  dose  *  week. 
However,  the  whole  plot  part  of  the  analysis  is 
erroneous.   Combining  the  whole  plot  part  of  the  split- 
plot  analysis  with  the  subplot  part  of  the  GLS  analysis 
produces  the  following  ANOVA  table: 


Table  7 

Source  r>f  Var  .  d  f 

Dose  4 

Error(Rat)  45 

Week  10 

Dose  *  Week  40 

Error(Week)  450 


ANOVA 


SS 


10295.72 

75668.30 

78626.58 

787.88 

4182.06 


MS 


2573.93 

1681.52 

7862.66 

19.70 

9.29 


-  E. 


1.53 

846.04 
2.12 


33 


The  conclusions  supported  by  the  GLS  analysis  in 
this  case  do  not  differ  from  the  conclusions  supported 
by  a  standard  split-plot  analysis.   In  either  case  there 
is  a  significant  dose  by  week  interaction  at  a  0.05 
level  of  significance.   However,  it  is  worth  noting  that 
the  week  sums  of  squares  dropped  substantially  from 
1783.51  to  846.04  . 

To  do  this  analysis  a  550  x  121  design  matrix  had 
to  be  produced.   A  laborious  task  to  say  the  least.   The 
final  proposed  method  eliminates  the  need  to  produce  the 
design  matrix  and  is  therefore  advocated.   SAS's 
generalized  linear  model  (GLM),  for  example,  will 
produce  the  design  matrix. 

GLM-GLS 

To  illustrate  GLM-GLS,  let  us  review  the  rat  data 
once  more.   As  with  GLS,  the  first  step  involves  running 
a  split-plot  analysis  of  the  data  using  a  procedure  such 
as  SAS's  PROC  GLM.   The  F-statistic  from  this  analysis 
will  be  correct  for  the  whole  plot  treatment,  dose  (see 
table  5  above ) . 

+ 
The  next  step  involves  converting  y    to  y    by 

^^^  i^k 

a  direct  transformation. 

The  final  step  involves  running  a  weighted  least 
squares  (WLS)  on  the  121  variables  of  the  split-plot 
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design  matrix  (X)  created,  for  example,  by  SAS's  GLM 

procedure  and  on  the  variable  of  the  transformed 

+ 
response  vector,  y  .   The  matrix  of  weights  to  apply  to 

WLS  is  of  the  form  W  =  I   0  G  ,  where, 

an 


4  0  0 
0  10 
0    0    1 


0  0  0 
0  0  0 
0    0    0 


0    0 
0    0 


0 
0 


0    10 
0    0    1 


11x11 


The  F-statistics  from  this  analysis  will  be  correct 
for  the  subplot  treatment,  week,  and  for  the  interaction 
term,  dose  *  week.  However,  the  whole  plot  part  of  the 
analysis  is  erroneous.  Combining  the  whole  plot  part  of 
the  split-plot  analysis  with  the  subplot  part  of  the  GLS 
analysis  produces  the  following  ANOVA  table: 


Table  8 


ANOVA 


Source  of  Var . 

df 

SS 

MS 

F 

Dose 

4 

10295.72 

2573.93 

1.53 

Error (Rat) 

45 

75668.30 

1681.52 

Week 

10 

78626.58 

7862.66 

846.04 

Dose  *  Week 

40 

787.88 

19.70 

2.12 

Error(Week) 

450 

4182.06 

9.29 
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Note:  there  is  no  difference  in  the  ANOVA  table 
presented  from  the  GLS  analysis  versus  the  table 
presented  from  the  GLM-GLS  analysis.   However,  the 
parameters  (0)  that  are  estimated  by  the  respective 
methods  differ  drastically. 

CONTRASTS 

SAS  uses  set-to-zero  restrictions  to  compute  the 
GLM  and  GLM-GLS  parameter  estimates.   Therefore,  the  GLS 
parameter  estimates  outputted  by  SAS  is  of  the  following 
form: 


CT  = 


T  -  T 

1  t 

T  -  T 

2  t 

T  -  T 

3  t 


T    -  T 
t-1    t 


To  illustrate  the  procedure  for  obtaining  GLS 
contrasts  from   the  GLM-GLS  parameter  estimates,  let  us 
take  one  last  look  at  the  rat  data.   Table  9  below  shows 
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the  parameter  estimates  for  both  GLS  and  GLM-GLS  for  the 
rat  data.   These  are  estimates  for  the  subplot 
treatment,  week. 

Table  9:   GLS  and  GLM-GLS  Parameter  Estimates 


4    11 


5    11 


T   -  T  =  -  36.18 
6    11 


T   -  T   =  -  29.06 
7     11 


GLS  GLM-GLS 


+  "   + 
T   -  T   =  -  66.20  T   -  T   =  -  30.788 

1  11  1     11 

T   -  T   =  -  60.72  T   -  T   =  -  25.308 

2  11  2     11 


+    '•      + 
T   -  T   =  -  54.80  T   -  T    =  -  22.676 

3     11  3     11 


T   -  T   =  -  49.14  T   -  T   =  -  20.568 


4     11 


+  "   + 
T   -  T   =  -  42.88  T   -  T   =  -  17.704 


5     11 


+  "   f 

T   -  T   =  -  14.760 
6     11 


+  ^   + 
T   -  T   =  -  11.660 
7     11 


+  "   + 
T^  -  T   -  -  22.52  T   -  T   =   -  9.392 


8     11 


8     11 
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GLS  and  GLM-GLS  Parameter  Estimates 


GLS 


T  -  T  =  -  15.32 
9    11 


GLM-GLS 


+  '^   + 
T   -  T   =   -  6.116 
9    11 


T    -  T   =  -  7.18 
10     11 


+   "   + 
T    -  T   =  -  2.296 
10     11 


0.00 


0.00 


The  desired  contrasts,  CT  ,  for  the  subplot 
treatment,  week,  can  be  obtained  by  premultlplying  the 


*+       -1 
GLM-GLS  estimate,  T  ,  by  CD 


For  example. 


the  contrast  T   -  T  can  be  obtained  by  multiplying; 
1    2 


-[ 


CT  =1-1   000000000 


-  66.21 

-  60.72 

-  54.80 

-  49.14 

-  42.88 

-  36.18 

-  29.06 

-  22.52 

-  15.32 

-  7.18 

0 
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-1 


=      CD 


T      =       -    5.49 


where. 


-1 


2.5  0  0  0  0  0 

1.5  10  0  0  0 

0.9  0.6  10  0  0 

0.54  0.36  0.6      10  0 

0.32  0.22  0.36    0.6  1  0 

0.19  0.13  0.22    0.36  0.6  1 

0.12  0.08  0.13    0.22  0.36  0.6 


0 
0 
0 
0 
0 
0 
1 


0 
0 
0 
0 
0 
0 
0 


0.07    0.05    0.08    0.13    0.22    0.36    0.6      1 


0 
0 
0 
0 
0 
0 
0 
0 


0.04    0.03    0.05    0.08    0.13    0.22    0.36    0.6      1 


0 
0 
0 
0 
0 
0 
0 
0 
0 


0.03    0.02    0.03    0.05    0.08    0.13    0.22    0.36    0.6      1 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 


0.02    0.01    0.02    0.03    0.05    0.08    0.13    0.22    0.36    0.6    1 


-  30.788 

-  25.308 

-  22.676 

-  20.568 

T   = 

-  17.704 

-  14.760 

-  11.660 

-   9.392 

-   6.116 

-   2.296 

0 
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-1  "  + 
Note:   CT  may  not  exactly  equal  CD   T  due  to  rounding 

error . 

In  summary,  three  techniques  were  presented  for 
analyzing  the  subplot  of  a  repeated  measures  experi- 
mental design  with  an  autoregressive  error  structure. 
The  first  technique,  proposed  by  Box  (1954),  involved 
adjusting  the  degrees  of  freedom  associated  with  the 
subplot  treatment,  whole  plot  treatment  by  subplot 
treatment  interaction,  and  the  subplot  error. 

The  second  technique,  GLS,  involved  transforming 
the  design  matrix  and  the  response  vector.   The 
transformed  response  vector  is  then  regressed  on  the 
transformed  design  matrix. 

The  third  and  final  technique,  GLM-GLS,  requires 
regressing  a  transformed  response  vector  on  the 
original  design  matrix. 

In  deciding  on  which  technique  to  use,  it  is 
helpful  to  note  that  GLS  provides  tests  that  are  more 
powerful  than  Box's  proposed  technique.   However,  the 
design  matrix  needs  to  be  constructed  to  do  a  GLS 
analysis.   This  is  time  consuming.   However,  one  can 
avoid  the  laborious  task  of  creating  the  design  matrix 
and  still  retain  the  nice  properties  of  GLS  by 
performing  A  GLM-GLS  analysis. 
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Several  statistical  computer  packages  exist  that 
have  a  generalized  linear  model  procedure  with  a 
weighted  least  squares  option.   Such  packages  can  be 
used  to  perform  a  modified  generalized  least  squares 
subplot  analysis  of  a  repeated  measures  experimental 
design  with  autoregressive  subplot  disturbances.   One 
need  only  regress  an  appropriately  transformed  response 
vector  on  the  repeated  measures  design  matrix  using  the 
weighted  least  squares  option.   The  weights  need  to  be 
specified . 

The  modified  generalized  least  squares  estimators 
retain  all  the  properties  of  the  standard  generalized 
least  squares  estimators.   However,  a  standard 
generalized  least  squares  analysis  requires  the 
specification  of  a  design  matrix.   Computer  packages 
with  the  generalized  linear  model-weighted  least  squares 
option  produce  the  design  matrix  at  a  considerable  time 
savings  to  the  analyst. 

Contrasts  can  be  constructed  for  the  modified 
generalized  least  squares  technique  by  taking  advantage 
of  the  linear  relationship  that  exist  between  the 
modified  generalized  least  squares  estimators  and  the 
standard  generalized  least  squares  estimators. 


